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The next generation of ground-based gravitational wave detectors may detect a few mergers of 
comparable-mass M ~ 100— IOOOA/q ( "intermediate-mass" , or IMBH) spinning black holes. Black 
hole spin is known to have a significant impact on the orbit, merger signal, and post-merger ringdown 
of any binary with non-negligible spin. In particular, the detection volume for spinning binaries 
depends significantly on the component black hole spins. We provide a fit to the single-detector and 
isotropic-network detection volume versus (total) mass and arbitrary spin for equal-mass binaries. 
Our analysis assumes matched filtering to all significant available waveform power (up to I = 6 
available for fitting, but only I < 4 significant) estimated by an array of 64 numerical simulations 
with component spins as large as S\,2/M 2 < 0.8. We provide a spin-dependent estimate of our 
uncertainty, up to Sx^jM 2 < 1. For the initial (advanced) LIGO detector, our fits are reliable for 
M e [100,500]A/ o (M € [100, 16OO]M ). In the online version of this article, we also provide fits 
assuming incomplete information, such as the neglect of higher-order harmonics. We briefly discuss 
how a strong selection bias towards aligned spins influences the interpretation of future gravitational 
wave detections of IMBH-IMBH mergers. 



I. INTRODUCTION 



Ground-based gravitational wave detectors like LIGO 
and Virgo are presently taking data at and beyond de- 
sign sensitivity [TJ [2]. Over the next several years as 
upgrades are performed, these detectors' sensitivity will 
increase substantially [3]. Advanced detectors are very 
likely to see many few-stellar-mass black hole binaries 
formed through isolated [U [5] and dynamical jBHH] pro- 
cesses. Additionally, advanced detectors could see the 
merger signature of two intermediate-mass black holes 
(each M £ [100, 10 3 ]), binaries which might be formed 
in dense globular clusters [£F. Unless astrophysical pro- 
cesses strongly suppress black hole spin, that spin will 
have a substantial effect on all the components of the 
signal to which these ground-based detectors are sensi- 
tive: the late-time inspiral, merger signal, and (through 
the final BH spin) ringdown. Though analytic approx- 
imations exist to describe the early-time (inspiral) and 
late-time (ringdown) behavior of a spinning BH-BH bi- 
nary, the merger signal must be obtained numerically, in 
principle for all possible mass (rti\,m-i) and spin (Si, S2) 
combinations. Of the types of black holes likely to be de- 
tected in the near future, intermediate-mass black holes 
have masses and spins such that their entire (short) de- 
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tectable waveform is dominated by the merger signal. As 
a result, very few candidate intermediate- mass merger 
waveforms are presently available, particularly for generic 
spins |10H13j . Conversely, the performance of real gravi- 
tational wave search pipelines is difficult to assess without 
extensive Monte Carlo simulations. All searches suffer 
from highly nongaussian noise and time-variable detec- 
tor performance; matched filter searches adopt a range of 
approximate waveforms, coincidence, and nongaussian- 
noise rejection strategies. However, given the compu- 
tational burden of each NR simulation and the limited 
selection currently available, extensive Monte Carlo stud- 
ies are not presently practical. Nonetheless, all numeri- 
cal simulations to date suggest that gravitational merger 
waveforms are surprisingly simple. For example, for 
aligned spin both the final spin [T4lfT7] and even merger 
waveforms [181 - 120"] have been accurately fit, for all pos- 
sible component masses (mi, 7712) and spin magnitudes 
I Si |, IS2I for moderate spin magnitude (< 0.9) and mass 
ratio (< 1/10), for the dominant mode of radiation. 



Given the simplicity of numerical merger waveforms 
and the few available simulations, in this paper we outline 
a simple method to analytically estimate the performance 
of present and future gravitational wave searches, extrap- 
olating from a small array of existing numerical simula- 
tions. Our method relies only on the raw numerical simu- 
lation output and detector response functions; we neither 
model the waveform itself nor limit to a few "dominant 
harmonics" associated with the early- or late-time binary 



orientation. 1 In Sec. Owe describe (i) how we estimate 
the detection volume of present and future gravitational 
wave searches for comparable mass binaries with known 
waveforms and (ii) how we extrapolate between them us- 
ing fits. In Sec. |HI| and Tab ic |Tl| we describe the set of 
NR simulations used. In Sec. |IV[ we compare our results 
on aligned spin to previously published estimates. Addi- 
tionally, we use our aligned spin results to emphasize how 
sensitive our predictions are to small systematic issues 
such as wave extraction radius and (to a lesser extent) 
numerical resolution. Then, in Sec. [V] we provide the 
coefficient functions needed for arbitrary spins with total 
and incomplete (l max — 2,3,4,...) waveform catalogs, 
along with our best estimates for parameter-dependent 
uncertainty in the detection volume. Finally, in Sec. |VI| 
we discuss how much (or little) spin influences searches 
for high-mass IMBH-IMBH mergers. 



II. GW SEARCHES FOR HIGH-MASS 
MERGERS: SELECTION BIASES 



The sensitivity of a network of gravitational wave de- 
tectors to a single class of randomly oriented source is 
often characterized in three ways: (A) via the maximum 
amplitude an optimally oriented but otherwise identi- 
cal binary produces; (B) using the angle-averaged signal 
power p* incident on a single detector [2T] from sources 
at a fixed distance; (C) via the expected detection rate for 
that class of source [4], adopting a distribution p(X) of 
sources described by parameters A (= component masses 
TOi,TO2; spins Si,S2; emission direction n in the frame 
of the binary; sky location and polarization angle N,ip; 
and distance r). For example, for a single interferometer 
with Gaussian noise, the signal-to-noise ratio p*(A) can 
be expressed as 

h det {t) = F + (N,-^)h+{h,t) + F x {N,-^j)h x (h,t) 
Mf)*b(f) 



(a\b) 



d.f~ 



S h (f) 



Pi = {hdet\hdet) 



(1) 

(2) 



where i*+,x are standard single-detector beampattcrn 
functions and where we use a subscript * when refer- 
ring to a single detector; compare to Appendix [C| For the 
first method, a peak signal to noise /5* imax , is particularly 
useful for aligned or nonspinning binaries dominated by 
I — \m\ — 2 emission, where the relative change in p ver- 
sus N,ip,ii is known analytically; see, e.g., Eq. (2) in 
[3]. In the second method, the source- and sky-location- 
averaged signal power incident on a single detector leads 
naturally to orientation-averaged power over a complex 



1 In other words, we neither fit to h of t or f nor do we pick a single 
"dominant mode" like /122, which relies on a preferred frame. 



wave amplitude 
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where p 2 is a technically convenient lower-dimensional 
average with clear physical meaning - the average over 
all source orientations of the network signal-to-noise re- 
covered by either (i) a pair of identical detectors, ori- 
ented at 45° to each other and with the source directly 
overhead, or equivalently (ii) a network of detectors with 
equal sensitivity to both polarizations in all directions 
[cf. Appendix [C] : 



dQ n 



(h + + ih x \h+ + ih x ) 



x / fixed h ' 



(5) 



Thus, the orientation-averaged power is technically con- 
venient since, if the emitted waveform is expressed as an 
expansion of the asymptotic complex waveform h + +ih x 
or curvature scalar ^4 into spin-weighted spheroidal har- 
monics 



h, 



ih y 
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then the angle-averaged signal power becomes a sum over 
inner products of the harmonic amplitude functions hi m 
or (with a different inner product) ^4,,i m given by 
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see, for example, Eqs. (7,8) in Reisswig et al. [19]. For 
clarity we have described both methods (A) and (B) as 
a characteristic SNR p for sources at a fixed distance, 
adopting a detection-strategy-neutral characterization. 
If one adopts a fiducial signal to noise ratio p c , such as 
a cutoff for single-detector SNR, then these amplitudes 
convert to physical distances; for instance, (B) implies 
an angle- averaged reach D* defined by the solution to 
p*(X, D) = p Cl or equivalently by 



P*W 



Pc 



A. (A) 



(9) 



where A = (mi, 7712, Si, S2 and suitable orbital phases) 
are the physical ("intrinsic") parameters of the binary 
(i.e., all parameters except %j},N,n and distance). 



Though well-defined, both methods only approximate 
the astrophysically-relevant sensitivity of gravitational 
wave detectors to spinning systems. For example, both 
the distribution and even optimal emission direction 
(-0, n) depend strongly on the direction of the total an- 
gular momentum J in band; and, therefore, on the 
masses and spins involved. While the peak amplitude 
could easily be tabulated and fit following the proce- 
dure described below, for astrophysical purposes (A) and 
(B) alone lose information about the beampattern shape 
function needed to construct (C), the astrophysically- 
relevant measure of sensitivity. To encapsulate all needed 
information about the beampattern shape, we introduce 
a beampattern function u>* for the ratio of single-detector 
SNR p to orientation-averaged single-detector SNR p*: 



P* = 



= p*(A)/p*(A) 



(10) 

(11) 



By construction, the orientation average of w% is always 
exactly unity ((iu«) = 1). In terms of this beampattern 
function and the previously defined angle-averaged range 
Z), the detection rate for sources in the nearby universe 
can be expressed as a sum over the rate per unit phys- 
ical volume dV = r 2 drdVL^ and per volume in binary 
parameters dX = dfl source d\: 



R 



D 



dN 



p>p, 
dN 



dtdVdX 



dVdX 



dtdV 
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source J p^ >p c 



dV. 



The last factors represent the detection volume averaged 
over source orientations. In the nearby universe, this 
orientation-averaged detection volume can be evaluated, 
reducing the detection rate to 



dN f ,,, „ 47T.= ,, 



(12) 



which shows that the final detection rate can be explained 
as a product of: (i) total event rate per unit volume 
jjdv'' ( n ) distribution of events in parameters; and a (iii) 
source-frame-averaged volume that characterizes the typ- 
ical reach to sources with parameters A, consisting of (iv) 
an orientation-averaged range D t , proportional to the 
band-limited SNR p* for a source at a fixed distance, 
and (v) a beaming correction factor w, , that depends on 
how the binary's polarized, beamed emission interacts 
with our network's polarization-dependent beampattern. 
In other words, the angle-averaged reach D* almost de- 
scribes the astrophysically relevant reach, modulo a weak 
correction factor 



u>*(A) 



r d£l n d£lNdil> 
J 47r(47r)(ir) 
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1/3 



(trt 



\2/3 



(13) 



which, like p, can be tabulated and fit versus all intrinsic 
parameters A. As we will see below, the beampattern 



correction function average u>* is necessarily almost al- 
ways unity, with ui* — 1 largest for aligned binaries and 
exactly zero for isotropic emission. In fact, given the un- 
certainties expected in our fit to /5(A), an excellent first 
approximation to u)* is unity. 



Fit versus spin 

Owing to the relative simplicity of numerical merger 
waveforms, scalar functionals of the waveforms like the 
final mass My and spin J 2 can be fit across the entire 
space of intrinsic parameters A 22J. Fits for the rem- 
nant BH's recoil kick and spin have been extensively 
explored in the literature [HI [15l [22H26] , Experience 
from fitting the final spins and signal-to-noise-ratio of 
aligned spinning binaries [Sec. IV suggests an accurate 
fit to p(Si, S2) requires at least cubic order in the com- 



ponents of Si,2 along the early-time orbital angular mo- 
mentum direction L (henceforth denoted z). 2 Though 
a generic symmetry-preserving expansion contains many 
components (roughly 6 3 /4 in a cubic-order expansion of 
a six-dimensional space with two Z2 symmetries, par- 
ity and black hole exchange 3 ), the physics of a process- 
ing merging binary strongly suppresses most terms. A 
truly generic symmetry preserving expansion of a scalar 
function allows scalar functions of Si, S2 to have a pre- 
ferred spin direction perpendicular to the total angular 
momentum, along the axis connecting the two holes at 
our simulations' starting time. Precession rapidly evolves 
in-plane spin components away from this preferred axis. 
We anticipate and test simulations confirm minimal de- 
pendence of the amplitude p on the relative orientation 
of spins to this preferred axis. 

For this reason, rather than use all spin components, 
as in Boyle and Kesden |22j . we describe our expansion 
in terms of three quantities: the preferred out of plane 
direction z ex J; the in-plane projection operator P per- 
pendicular to z; and a pair even- and odd- exchange- 
symmetric reduced spins: 

X± = (mid ± m 2 a 2 )/Af [S oc -%_] (14) 

where a& = Sk/m 2 . An exchange-symmetric expansion 
of equal-mass binaries must have even powers of X- • [For 
unequal masses, terms proportional to X- are possible 
when prefixed by asymmetric mass terms. A generic ex- 
pansion will introduce at linear spin order (5m/M)(x~ ' 
z), at quadratic spin order (6m/M)Px- -Px+: el cetera. 



To avoid ambiguity. I adopt the initial orbital ang ular momentum 
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L at r = 6.2 as a reference direction; see Section '. 
Following 1221 a generic scalar expansion would have 1,2,11, and 
23 parameters at zeroth, first, second, and cubic orders, respec- 
tively. Not including the trivial zeroth order term, our proposed 
expansion has 4 fewer parameters at quadratic order and 12 fewer 
parameters at cubic order. 



The exchange-symmetric coefficients provided below can 
depend on only even powers of (5m/ M).] Working to 
cubic order, we anticipate p has the form 

p(m 1 ,m 2 ,S 1 ,S 2 ) = p (m 1 ,m 2 )[l (15) 

+ X 1 ( X+ ■ z) + X 2 ( X+ ■ if + x s (x+ ■ z) 2 

+ X 02 (P X+ ) 2 + A 20 (x„ ■ zf + Aai{PX-f 

+ Bi2 o(x+ ■ z)(P X+ ) 2 

+ B 102Q (x+ ■ z)(x- ■ zf 

+ B 1002 ( X+ • z)(P X -) 2 + 0(X 4 ) • • •] 
uj*(mi,m 2 ,Si,S 2 ) = uj ,,(mi,m 2 )[l (16) 

+ ^ix+ -z + Z 2 (x + -zf... 
+ C 20 (x- ■ z) 2 + C 02 {PX-) 2 + ■ ■ ■ 
+ V 1200 ( X + ■ z){PX-) 2 + ■ ■ ■ 
= 1.095[l + ^ix+-z + ---] 

where p (mi,rri2) is provided by a comparable calcula- 
tion for nonspinning binaries and where the value for u) * 
for Si = S2 = shown can be estimated using only the 
dominant I = \m\ = 2 aligned emission beampattern. 
In most of the text we will refer explicitly to the coeffi- 
cient functions described above. However, when referring 
to this expansion in its entirety, we adopt the shorthand 
abstract notation y a for its coefficient functions (of M) 
and ip a for basis functions (of x±) : 



y a in Eq. (15), adopting uniform uncertainties in p k for 



P/Po = ^2 ip a y a 



(17) 



For equal-mass binaries mi = m 2 = M/2, we deter- 
mine the coefficient functions X\,2,3 (M), . . . suitable to a 
specific gravitational wave detector noise power spectrum 
Sh as follows. 5 We pick a total mass M. For all numerical 
simulation k in Table [TTl corresponding to spin combina- 
tions <Si,fc, S^fc, we extract time domain spin-weighted 
harmonics of the Weyl scalar ^A.im,k{t) for I < 6. We 
construct each p k using Eq. pi). Likewise, we construct 
w by (i) reconstructing h along a large number of ran- 
domly chosen orientations (N,i/^,n), using the multipole 
coefficients ^A,lm] (ii) calculating w^; then (iii) averaging 
over all the random samples. Finally, excepting only a 
few simulations chosen as blind tests of our fit [Sec. M, 
we perform a simple least-squares fit 6 for the coefficients 



4 To an excellent approximation the beampattern factor w agrees 
with the average associated with I = \m\ = 2 emission and can be 
calculated by angle-averaging the analytically-known form of to 
for that case; see O'Shaughnessy et al. [4], noting w there differs 
by a constant factor from our definition. 

5 We adopt an the initial LIGO sensitivity from 1271 and an ad- 
vanced detector noise power spectrum from [3]- 

6 Though we prefer a global fit to all parameters, the individ- 
ual coefficient functions themselves can be (nearly) isolated us- 
ing suitable-symmetry one-parameter subfamilies of simulations. 
For example, the coefficient A 2 o can be determined from simu- 
lations with antialigncd spins Si = — S 2 = |5|z; the coefficient 



all fc. Since w varies little from unity and since our esti- 
mate is subject to significant Poisson sampling error, we 
retain only leading-order dependence with spin (Z\). 
Fit errors: Numerical, systematic, and truncation errors 
limit our ability to determine these coefficients precisely. 
High-order or symmetry-suppressed coefficients like X 3 
and A20 are particularly sensitive to small numerical er- 
rors. Furthermore, the recovered coefficients y a (M) have 
highly correlated uncertainties E Qj g(M), where E a ^ is 
the least-squares estimate of the parameter covariancc 
matrix. Finally, phenomenologically speaking the inter- 
esting uncertainty is how much our fit F to our data, say 
for p, 

Fk = Pk/Po , 

differs from truth. Though highly dependent on the pa- 
rameter distribution of simulations used in the fit, one 
simple estimate for overall error is the residual rms error 
between the data and our least-squares fit: 



SF 2 = J2(F(X k )-F k ) 2 /N. 



A far more stable and spin-magnitude-dependent repre- 
sentation of fit uncertainty is the expected L 2 error F(a) 
in the fit for spins with magnitude Si/m 2 , S2/1TI2 < a: 



T 2 (a) = ((p-(p)) 2 ) m /p 2 



^ E ab (lp a 1pb) 



,s a + s b 



(18) 



where we use integer exponents s a to describe the order 
of the basis functions ( / a (xSi, XS2) = x Sa ip a (Si, S2)); 
where S Q ^ is the least-squares covariance matrix of 
fit parameters; and where the coefficients are easily- 
tabulated moments of the coefficient matrix over all 
spins, provided in Table llj 

{VaWau = / ,. /on 3 Va(a u a 2 )^ b (a 1 , a 2 ) .(19) 

Roughly speaking, .F(a) estimates the relative error in 
p when applying our fit to a generic pair of spins with 
typical magnitude a. 

Astrophysical tolerance: Astrophysically speaking, the 
orientation-averaged range relates the number of sources 
observed (or their absence) to the implied source event 
rate (or upper bound). Even adopting the most opti- 
mistic assumptions, no more than O(10) intermediate- 
mass mergers are expected to be seen by advanced 



61020 from simulations with aligned but unequal spins, given 
knowledge of equal-spin coefficients; the coefficient .4o2 can be 
determined from the "B-series" |28| . where both spins are an- 
tialigncd (x+ = 0) and tilted at a range of angles 8. For brevity, 
we only discuss a global fit, rather than fits to individual simu- 
lation subfamilies. 



X\ X-2 Xz X()2 A20 A02 B1200 S1020 ^1002 

0.168 0. 0.065 0. 0. 0. 0.047 0. 0.047 

X 1 0. 0.065 0. 0.047 0. 0.047 0. 0. 0. 

Xi 0.065 0. 0.032 0. 0. 0. 0.015 0. 0.015 

X 3 0. 0.047 0. 0.175 0.046 0.064 0. 0. 0. 

X 02 0. 0. 0. 0.046 0.063 0.047 0. 0. 0. 

A20 0. 0.047 0. 0.064 0.047 0.176 0. 0. 0. 

A02 0.047 0. 0.015 0. 0. 0. 0.021 0. 0. 

£1200 0. 0. 0. 0. 0. 0. 0. 0. 0. 

£1020 0.047 0. 0.015 0. 0. 0. 0. 0. 0.021 

TABLE I: Coefficient matrix (ip a '>pb) described in Eqs. 
( T8|l9) . 
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FIG. 1: Top panel: Expected la error T in an unconstrained 
fit to F — p/p , for spin magnitudes |ai | , 1 0-2 1 limited to below 
unity (blue; largest); below 0.8 (red); and below 0.6,0.4 and 
0.2 (similar curves near bottom). Also shown (points) are the 
least-squares errors SF between our fit F and the simulated 
NR data. Bottom panel: As above, but using a model where 
only Ai : 2,3 and A02 are nonzero. 



was used in previous binary black hole (BBH) studies 
[28H36] . The grid structure for each run consisted of 
10 levels of refinement provided by CARPET [57], a mesh 
refinement package for CACTUS [38 . Sixth-order spatial 
finite differencing was used with the BSSN equations im- 
plemented with Kranc [39 . The outer boundaries are 
located at 317M. Each simulation was performed with a 
resolution of M/77 on the finest refinement level, with 
each successive level's resolution decreased by a factor of 
2. All BBH simulations have two equal-mass black holes 
(BHs) with total mass M = mi + m 2 = 1.0 initiated on 
the x-axis, with initial separation as in Table |TT| Table ITT] 
also lists the initial spin configuration and the length of 
the simulation and waveform. 

Because we attempt to fit to small changes in the am- 
plitude versus spins, we carefully estimate the effects of 
possible numerical artifacts, such as waveform extrac- 
tion radius and simulation resolution. For example, in 
a few cases waveforms were generated at alternate res- 
olutions; convergence consistent with our fourth order 
code is found. 

Our best fit coefficient functions and their (fitting) un- 
certainties are provided in Figures [3] and [H] These fits 
reproduce data from our highest-resolution simulations, 
extrapolated to infinite radius. 

As seen in Table [TT] we employ two choices for the 
initial separation: r = 10, 6.2M. The spins and orbit 
of binaries started at r = 10M precess, evolving into a 
slightly different configuration by r = 6.2M. The spin 
and L configuration at early times is not identical to 
(but can be reconstructed from) our simulations' start- 
ing point. Similar precession effects have been included 
in fits to merger recoil kicks, to correctly reconstruct the 
kick direction as a function of spins at very early times 
p~5l IT7] . These 9 spin-misaligned systems are dominated 
by their aligned spins and thus precess only very slightly 
between r = 10,6.2 (e.g., z ■ L > 0.96). However, to 
avoid systematic errors caused by different starting radii, 
for these simulations we extract the (coordinate) spins 
and L at r — 6.2. In fact, our overall answer changes 
little, independent of whether r = 6.2 or r = 10 is used. 



IV. COMPARISON WITH PREVIOUS 
RESULTS: ALIGNED SPIN 



ground-based detectors [3], at best allowing the source 
event rate be determined to O(30%) at ltr confidence. To 
achieve this level of accuracy, we require only O(10%) ac- 
curacy in Dw. Even for maximally rotating black holes, 
our fits should be at least that accurate; see Fig[TJ 



III. NR SIMULATIONS 

Table [TTJ lists the waveforms used in the present work. 
These waveforms were produced with MayaKranc, which 



For equal-mass binaries with spins aligned with the 
orbital angular momentum (Px+ — PX- = 0), many 
previous studies have provided numerical waveforms [401 - 
132] i detailed tabulations of the gravitational wave ampli- 
tude p available to different gravitational wave detectors 
|19) . mismatch-based estimates of waveform complexity 
[4"31 |4"4"] . and even phenomenological fits to the wave- 
forms themselves [15] P4"5l447| . To adopt a fiducial ref- 
erence which directly provides comparable information, 
we compare our fits to generic spins to the spin-aligned 
results provided by Reisswig et al. [IT?] . 

Equal-mass spin-aligned systems preserve a symme- 
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TABLE II: The set of equal-mass merger simulations used in this paper. As described in Section |Vj the first 2 simulations 
have randomly chosen spin orientations; were not used in any fit; and provided a blind test of our fitting procedure. In this 
table, initial conditions are specified by the first six columns, which provide the component spins Sk/M 2 of each black hole, 
along with the initial separation r start- Two columns (T,T wave ) provide the duration of the simulation in its entirety and of 
the resolved, conversing portion of the I = m — 2 waveform, respectively. Finally, hmin shows the smallest resolution used in 
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FIG. 2: Angle averaged SNR p (note: not p* = p/s/S) ex- 
pected for the present LIGO detectors, extracted from our 
numerical simulations (solid lines) and from our fits (dotted 
lines indicate our preferred estimate from an unconstrained 
fit to all simulations; shaded interval indicates our estimate 
of fitting error .F(0.8)po [Eq. [18]). The colored solid lines 
correspond to equal-mass binaries with a\ — a 2 = 0.8 (blue, 
red), placed at r = 100 Mpc, based on matched filtering to 
a single-detector datastream containing only the harmonics 
(i) / = 2 (blue, solid); (ii) I < 4 (blue, dotted). For com- 
parison, the thin black (dotted) line indicates our estimate of 
p(M) using I < 2 (I < 4) for a binary of nonspinmng black 
holes at the same distance. Compare with Fig. 6 of Reis- 
swig et al. [T§] and our Figure |4l noting (i) p — p*\/5 and 
p* = p*max/(5/2) for the / = 2 subspace (blue) and (ii) for 
the I — 4 contribution they show the peak magnitude p #! max 
(depending linearly on the amplitude due to the I = 4 sub- 
space) while we show the angle-averaged magnitude p (de- 
pending quadratically on this amplitude and therefore highly 
suppressed). Both papers adopt a comparable initial LIGO 
noise curve. At those M with a large detection volume, the 
SNR including higher harmonics is nearly indistinguishable 
from the SNR from I — 2 alone. Conversely, for sufficiently 
high masses (e.g., M > 600Mq for a — 0), the I = 4 mode 
dominates the angle- averaged power p; see, for example, the 
two curves corresponding to a = 0. However, though I = 4 
does dominate at the highest masses, systematic errors as- 
sociated with extrapolating the extraction radius to infinity 
can lead to apparent contradictions. Here, for a = 0.8 the 
extrapolated "total" SNR from all modes I < 4 is slightly less 
than the corresponding extrapolated p(M) curve including 
only 1 = 2. 



FIG. 3: Fit coefficient functions #1,2,3 (blue, red, black, 
respectively), A20 (green), and ,61020 (orange), along with 
1<t least-squares fitting uncertainties %/^fcfc m each param- 
eter. [Error bars are suppressed for the ill-constrained cubic- 
order term ,61020 (orange). This cubic term's errors large and 
strongly correlated with others; also, this term's preferred fit 
evolves significantly with extraction radius, even though each 
slice is consistent with 61020 ~ 0.] Solid lines show param- 
eters of our general multiparameter fit to all available sim- 
ulations; dotted lines correspond to fits to only aligned-spin 
data 5*1,2 oc z using a special restricted model only X\o, are 
nonzero. 



eraged amplitude p* are related by 



P*max r. P* 



1=2 only 



(20) 



The I > 2 modes in general and the / = 4 modes in par- 
ticular violate this relation: the signal power p* tm ax seen 
viewing along the symmetry axis increases linearly with 
I > 2 amplitude; however, the orientation-averaged power 
p increases quadratically in the higher- harmonic ampli- 
tude (here, 1 — 4). Keeping this distinction in mind, 
[Fig. [2] agree with Reisswig et al. [TS] [their Fig. 6], 
who find that for lower-mass mergers M < 5OOM ,1 = 4 
modes are typically a 5—10% correction in p max (depend- 
ing on spin) to the individual signal-to-noise ratios p max , 
but a much smaller correction to p and the ratio p/p - 
Moreover, as in all previous studies, aligned-spin wave- 
forms appear to depend only extremely weakly if at all 
on the antisymmetric spin combination (x_) |44j . mostly 
through power-suppressed higher harmonics, 8 suggesting 



try direction throughout their inspiral, insuring the I = 
\m\ = 2 modes dominate above higher (even) I orders. 
Limiting attention to the I = 2 subspace, the peak am- 
plitude p* t max and source- and detector- orientation av- 



Rcisswig et al. 19 provide a comparable expression, p t ,avg = 
\/5j2p*,max, relating the peak SNR p*, maa; to the source- 
orientation-averaged SNR p* } avg- Their paper adopts a similar 
notation, without a single-detector * subscript. 
Previous studies have demonstrated that the 2,2 mode of 
spin-antialigned systems resembles the nonspinning waveform. 
When higher harmonics are included, the nonspinning and spin- 
antialigned waveforms have appreciable differences, measured by 
their relative mismatch. However, these higher modes carry com- 
paratively little power: see, for example, Figure 2 of Shoemaker 



-^20 = S1020 ~ 0. Additionally, as previous studies have 
shown, the orientation-averaged amplitude p increases 
monotonically with x+ ' z - in terms of our expansion, 
its parameters must satisfy d x .±F > 0, or 



Xi+2X 2 ( X+ -z) + 3X a ( X+ -Z) 2 >0. 



(21) 



Figure [3] provides our best estimates for the coefficient 
functions for p relevant to aligned spins, both employing 
all data (solid) and only I = 2 aligned data (dashed); 
all satisfy this property. Though our general fit allows 
for small nonzero -420>£>i020> X 2 , our equal-mass aligned 
spin data is well fit assuming these three parameters are 
exactly zero and the aligned-spin model consists of only 
Xi t s (dotted lines in Figure [31 for comparison, the fitting 
error SF of this restricted model to aligned-only data 
is comparable to the fitting error shown in Figure [T] for 
a generic model to all data). Likewise, as our global 
fit error is often comparable to the change in p due to 
higher harmonics, particularly beyond I = 4 for M < 
5OOM0, these harmonics can to a first approximation 
be neglected. Finally, for the mass range where I = 2 
emission dominates, the beampattern is extremely well 
approximated by the fiducial nonspinning values. 

For a fixed bandpass detector, high mass favors higher 
harmonics, particularly when those harmonics are emit- 
ted from inspiral while the I — 2 modes arise from the 
exponentially suppressed ringdown phase. For aligned 
binaries above M ~ 500M Q , though the overall range 
drops, the I = 4 harmonics provide an increasingly sig- 
nificant fraction of detectable power, even in the absence 
of spin; see Figure |4] As a result, based on the raw simu- 
lation data the beampattern correction for initial LIGO 
detectors evolves from w* sw 1.095 at low mass (appro- 
priate to I = \m\ = 2) down to a relatively isotropic 
w* « 1.04 at M ~ 600 Mq during the transition be- 
tween I = 2 and I = 4 beampatterns before eventually 
rising back at higher masses (i.e., approaching the value 
1.085 appropriate to I — \m\ — 4; see Table III and Ap- 
pcndixIC]). In principle, the process continues with higher 
harmonics; in practice, however, the detection volume is 
sufficiently small that their contribution to detection is 
negligible. 

Though we estimate p as the nonspinning value p 
times a small "correction" F, at high mass and spin the 
raw numerical data requires extremely large F. Compar- 
ing Figure[2]with Figurefll when the ratio of nonspinning 
to spinning p is large, our fit breaks down. 



et al. [44) . keeping in mind modes contribute power to p/p in 
quadrature. Therefore, though p depends slightly on the anti- 
symmetric spin X- 1 to a leading-order approximation their effect 
can be neglected. 
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FIG. 4: The relative significance of signal power due to dif- 
ferent harmonics, shown for a nonspinning equal-mass binary 
seen by initial LIGO (top panel) and advanced LIGO (bottom 
panel), shown at an extraction radius r — 75. The distances 
between the curves indicate p 2 Ll the signal power due to all p 
harmonics. At each mass, p\ is normalized to the contribu- 
tion from all Z < 4 modes (p< 4 ) For advanced detectors, the 
quadrupole mode strongly dominates for most masses. Ad- 
vanced detectors have much better low-frequency frequency 
sensitivity and particularly a much less steep low-frequency 
limit. On the contrary, low-mass detectors have a steep low- 
frequency boundary, strongly filtering out the I = 2 mode 
once the final mass and spin lead to ringdown out of band. 



Convergence I: Extraction Radius 

Our low-resolution large-radius grid zones do not re- 
tain enough information to permit adequately accurate 
waveform extraction for all harmonics. Therefore, unlike 
Reisswig et al. 19J who extract at r = 160M and partic- 
ularly unlike extraction at J + |5S], we extract relatively 
close to the binary, at coordinate radii r = 40, 50, 60 
(and, when available, at r = 75). As seen in Figure [5l fi- 
nite extraction radius effects can compete with the small 
spin-dependent changes in p. 

For this reason, rather than adopt a single preferred 
extraction radius, we first extrapolate p{M) to infinity, 
then fit to the extrapolated p(M) data. Owing to insta- 
bilities in the extrapolation at high mass to p calculated 
using initial detector noise spectrum Sh, we do not trust 
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FIG. 5: For an equal-mass aligned-spin binary with ai — 
a-2 = 0.8, estimates for p (top panel) and p/p (bottom panel) 
seen by a single initial LIGO detector at D = 100 Mpc versus 
total binary mass M, based on waveforms extracted from our 
numerical merger simulations at coordinate extraction radii 
r = 40,50,60,75 (blue, red, yellow, green), and (thick solid) 
linearly extrapolated in 1/r to r — > oo. For high- mass mergers 
M > 500 Mq, the aligned-spin range is much larger than the 
nonspinning range (p/po > 2; compare to Figure |2l. We 
anticipate both our extrapolation to large radius and our fit 
to perform poorly when the ratio p/p is large. 
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FIG. 6: For an equal-mass binary with M = 3OOM0 and 
generic spins, estimates for the fitting parameters Afi.,2,3 and 
61020 to p/po, based on numerical waveforms extracted on 
spheres of radius r = 40, 50, 60, shown versus 1/r. Colors and 
symbols are associated as in Figure [3] This figure provides 
neither point error estimates nor extremely strong correlations 
between recovered parameters. Also shown at 1/r = are 
the best-fit parameters to p/po data that has been linearly 
extrapolated to r -> 00 based on these radial slices as well 
as to the often- available r = 75 slice. Finally, the light gray 
region indicates the range of extraction radii associated with 
our two lowest-resolution zones (i.e., between r = 80 and 
317); to avoid introducing systematic error, we extrapolate 
over radii associated with a single grid zone. 



radius is comparable to the semianalytic error estimate 
shown in Fig. IT] 



V. SELECTION BIASES FOR GENERIC, 
COMPARABLE-MASS MERGERS 



our fit for extremely high masses M > 5OOM . On the 
contrary, for advanced detectors which lack such a steep 
low-frequency cutoff, our fitting procedure works well to 
proportionally higher masses; see Appendix |A| 

Figure [6] also demonstrates that specific coefficient val- 
ues are sensitive to extrapolation, particularly the high- 
order (cubic) coefficients. This trend versus extraction 
radius demonstrates a key issue associated with any phe- 
nomenological fit: while we can present best-fit parame- 
ter estimates, these estimates are subject to strong cor- 
relations. 

Despite uncertain and highly correlated fitting parame- 
ters, the best fit functions extracted from each extraction 
radius and from the radially extrapolated data are both 
extremely consistent with one another and with decreas- 
ing differences as r increases. Specifically, performing 
Monte Carlo estimates of the 1? difference between the 
finite-radius fits and the fit to extrapolated data, we find 
the set of L 2 functional differences \\F r — F^W (i) is con- 
sistent with being proportional to 1/r and (ii) at each 



Even well before the epochs considered here, bina- 
ries with generic spins precess, breaking symmetry and 
distributing power among other I = 2 modes. The 
asymmetric merger process and ringdown favors excit- 
ing otherwise-suppressed higher harmonics (e.g., / = 
3), as well as the I = 4 modes present even without 
spin. If the merger is in band, corresponding to masses 
M <G [50, 500] M Q , higher harmonics contribute an in- 
creasing proportion of overall SNR, particularly at high 
spin; see, e.g., Fig. 2 in Shoemaker et al. [U]. As in 
the nonspinning case, if late stages of merger are de- 
tected, corresponding to masses M > 500M Q for initial 
and M > 1000M Q for advanced detectors, the beam- 
pattern can be significantly asymmetric. Unfortunately, 
though substantial beampattern asymmetries begin to 
occur at high mass and spin (w*(M) varying), our fit to 
the angle-averaged power (p/p ) is not accurate enough 
at these high masses to justify a detailed analysis. At 
lower mass, spin precession only weakly anisotropizes the 
beam. Comparing to our symmetry expansion [Eq. ( |15[ )], 
any (linear) leading-order spin-dependent corrections de- 



10 



pend on only aligned spin components; the linear term 
is without loss of generality determined by simulations 
of aligned merging binaries. As the beampatterns from 
aligned merging binaries are dominated by I = \m\ = 2 
emission, the beampattern correction factor for generic 
misaligned binaries u)* changes little from its nonspin- 
ning value. Though the beampattern changes shape sub- 
stantially due to precession, its effective volume is nearly 
unchanged. 

On the other hand, the average amplitude p varies sub- 
stantially with spin magnitude and orientation. Well- 
known results for aligned spins illustrate the relative im- 
pact that spin orientations can produce. As a function of 
spin orientations relative to L, the largest p occur with 
spins that are aligned with the orbit (%_ = 0; Si 2 oc z); 
the smallest p occurs with spins that are antialigned with 
the orbit (x- = 0; Si ,2 oc — z); and p is essentially un- 
changed if the spins are mutually antialigned (x+ = 0. 
both for Si 2 oc z and generally for all spin orientations 
X-)- In fact, for these these intuitively obvious condi- 



tions to hold, the expansion parameters in Eq. ( 15 1 must 
satisfy the following conditions: 



-420 

3 



A 



02 



*i±ix + i*2+ 7 ix+r*3 > 





£1200 |x+ 1' 



±|X+I#0: 



all of which our best-fit coefficients satisfy. As noted pre- 
viously, because we included cubic-order terms, our fit 
for p is also consistent with p being a positive and mono- 
tonically increasing function of equal, aligned spins; see 
Eq. (21). More generally, for arbitrary spin orientations 
our fit is positive-definite for at least M < 500M© and 
| <xx,2 1 < 0.8, as well as for selected spin configurations 
at higher mass. As in the nonspinning case, however, 
our fit breaks down for binaries with large aligned spins 
and high mass, as in these regions the correction factor 
F — p/p must be nearly (for spins mostly antialigned 
with the orbit) or much larger than unity (for component 
spins mostly aligned with the orbit). 

provide our best-fit coefficients to the 



Figures [3] and [8 



expansion of Eq. (15). These preferred values reproduce 
our best simulation resolutions, extrapolated to r-> oc, 
including all available harmonics. 9 Though the symme- 
try expansion of Eq. ( 15 ) permits more generic behav- 



ior with spin, our numerical results are well-fit with a 
far more restrictive form where only #1,2,3 and #02 are 
nonzero. 10 As shown by the bottom panel in Figure 111 



9 For completeness, in the online version of this article we provide 
data for and fits to the cases where only I < 2, 3, 4, 5 harmonics 
were used. Except for low spin, fitting errors are comparable to 
or larger than the contribution of higher harmonics. 
10 In fact, roughly speaking this restricted fits' coefficients are 
themselves related by a nearly mass-independent proportional- 
ity: X 1 ss X 3 « 0(2 - A)X 2 « 0(6 - 8)#02, both for initial 
(M S [1OO,5OO]M ) and advanced (M e [200, 1300] M ) LIGO 
detectors. 




FIG. 7: For a M — 200M Q binary, angle- averaged SNR ratio 
F — p/po (points); our generic, unconstrained cubic-order fit 
F to all simulations (solid line); and a fit where only #1,2,3 
and #02 can be nonzero (dashed lines). The data and fits 
are shown for spin magnitudes a — 0.4, 0.6, 0.8 (green, blue, 
red, respectively) and spins Si in the orbital plane (along the 
separation vector) and S2 tilted by angle 9 away from z. 
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FIG. 8: Fit coefficient functions #02, Am, .81200,1002 (black, 
red, blue, purple) and estimated least-squares fitting uncer- 
tainties Sfefe for the coefficients relevant to misaligned spin; 
compare with Figure 13] As in that figure, the large fitting 
errors (> 0.25) on cubic-order terms (,61200, S1002) are not 
shown. Solid lines show a general fit including all harmonics 
and spins. The dotted black line corresponds to the value of 
#02 assuming a model where only it and #1,2,3 are nonzero. 



this fit performs well averaged over all simulations be- 
cause our simulations mostly have 1 01,2 1 < 0-6 [Table [TT| , 
where cubic order corrections are still small. At the same 
time, the fit has more than enough parameters to explain 
the predominantly linear- and quadratic-order variation 
in p versus generic spin orientations; see for example Fig- 
ure [7) 

Blind test The first two simulations in Table [Til with 
random spins with |oi| = 1 0.2 1 = 0.6, were reserved as 
a blind test. As illustrated by Figure [9j for low masses 
M < 5OO-M our fit correctly recovers p/p to within 
roughly the la relative error J 7 , shown as the shaded 
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FIG. 9: A plot of F — p/po for the "blind test" spin configu- 
rations. Dots show the raw results from our numerical simu- 
lation; the curves and shaded region indicate results from our 
preferred fit. 



region on this plot. 



VI. ASTROPHYSICS 



Using the fits provided in this paper, we can use Eq. 



( 12 1 to determine the relative likelihood of observing dif- 



ferent spin magnitudes and alignments, given a progeni- 
tor parameter distribution p(X). As expected from stud- 
ies of aligned spins, ground-based gravitational wave de- 
tectors will be biased towards the detection of aligned 
and (to a lesser extent) large spins. For example, as seen 
in Figure 11 in the special case M — 200M Q , popula- 



tions of randomly oriented spins with identical magnitude 
l a i I = 1 02 1 = a arc much less likely to be seen than a com- 
parable population of spin-orbit aligned spins, and con- 
versely for spin-orbit antialigned spins a\ = a 2 = — az. 
However, populations of binaries with random spins (dot- 
ted line in Figure 11 ) have comparable average detection 



volumes as nonspinning binaries, with detection volume 
increasing only O(50%) if spin magnitudes as large as 
0.8 are allowed, only marginally larger than fitting error 
in the detection volume. In other words, in a popula- 
tion of random spins, any pair of spin magnitudes are 
roughly likely to be observed, with only a slight bias to- 
wards larger spin magnitudes. Despite the much larger 



detection volume for aligned spins, a priori perfectly 
aligned and large spins should be rare, suppressing the 
bias towards high spin magnitude. For similar reasons, in 
any population of random spin orientations, large aligned 
spins occur rarely enough and the bias towards large spin 
is small enough (typically less than a factor 5) that the as- 
sociated detected population will not be overabundant in 
tightly aligned spins unless the progenitor population is. 
For the purposes of illustration, suppose we assume any 
spin pair with both spins S\$ within 7r/4 of alignment 
with z (i.e., S\ : 2 ■ z > l/v2) will be significantly ampli- 
fied, independent of spin magnitude. Even in this unrc- 
alistically optimistic case, only [(1 — cos7r/4)/2] 2 ~ 2% 
of all random progenitor spin orientations could be am- 
plified. 

All the properties described above and exhibited in 
Figure [TT] can be understood analytically. For simplicity, 
let us adopt a fit with only ^1,2,3 and Xq 2 nonzero (dot- 
ted lines in Figures 3p). The detection volume can be 
approximated to quadratic order as 



K(S!,S 2 ) ex ^[1+3Ai( Xh 

+ 3X a2 (P X+ ) 2 . 



z) + 3(* 2 

o(x 3 )] 



•*a)(x H 



(22) 



Averaging the detection volume over all spin orientations 
(or any symmetric volumes in Si, S 2 ) eliminates terms of 
odd order, leaving only quadratic-order dependence on 



spin: 



<nSi,S 2 )) n ex ^[l + 3(*i 2 + * 2 )(( X+ -z) 2 ) 

+ 3X 02 ((P X+ ) 2 ) n + O( X ) 4 ]. 



(23) 



This expression reproduces the corresponding (blue dot- 
ted) curves in Fig. 11 Strictly speaking, our simu- 
lations and fit apply only to equal mass ratio. How- 
ever, given symmetry considerations, the leading-order 
spin- and £m/M-dependent corrections to the detection 
volume must be quadratic in the asymmetric mass ra- 
tio times a quadratic function of spins. Our estimate 
therefore applies to many comparable-mass IMBH-IMBH 
mergers as well. 

Some BH-BH binary populations have an intrinsic bias 
towards aligned spin, such as the products of binary evo- 
lution of extremely low-metallicity binary black holes. In 
these cases, spin-orbit misalignment encodes otherwise 
inaccessible information about the strength of supernova 
kicks on those very massive black holes. 12 For an intrin- 
sically aligned population, gravitational wave detectors 
are far more likely to observe large, tightly aligned spins. 



Figure 10 illustrates this effect using contours of the de- 
tection volume versus X + -z = (Si z +S2, z )/2 (the aligned 



For spins with identical spin magnitude |ai2| = a, the angle- 
averaged coefficients of this expansion are provided by Table IT1 
12 Very massive BHs should have supernova kicks strongly sup- 
pressed by fallback, or even prompt collapse. The complete ab- 
sence of spin-orbit misalignment would confirm this model. 
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FIG. 10: Contours of single-detector initial LIGO detection 
volume to an equal-mass m\ — mi — 200Mq binary versus 
the symmetric spin x+, adopting a simple fit with only #1,2,3 
and A02 nonzero. Contours are shown for V — 2 n V for n = 
—4, —3 ... 3, increasing from left to right on the graph, where 
V is the detection volume for zero spin. The shaded regions 
indicate physically unattainable %+ values (i.e., which must 
violate lai 2I < 1). 



component of spin) and {Px+) 2 (the perpendicular part 
of the symmetrized spin), adopting a fit with only ^1,2,3 
and A20 nonzero. 



Below roughly 40 Hz, the initial LIGO detector's de- 
sign sensitivity to merger waves (0(f~ 7 ^ 3 /Sh)) decreases 
rapidly. Some mergers yield a final-state black hole whose 
ringdown frequency is comparable to this frequency. As 
a result, small changes in mass or spin lead to large 
changes in the volume to which these black holes can 
be seen. In particular, for black holes of fixed mass, 
this low-frequency cutoff produces a noticable bias to- 
wards large, aligned spins. Unfortunately, our fits to p 
degrade precisely where this selection bias becomes im- 
portant. However, as seen in the bottom panel of Figure 
[TT] the bias towards large spins is already apparent at 
M = 500-Mq. For example, at a — 0.8, the blue solid 
line (randomly oriented spins) is twice as large as the 
nonspinning value, in contrast to only tens of percent 
higher at lower mass. At this and higher masses, the 
distance to which a nonspinning binary is visible can be 
substantially greater, up to a factor of 6 or more near 
AI = 600M Q ; see Fig. [5J This substantial increase in 
detection volume oc (p/p ) 3 ~ O(10 — 100) can partially 
compensate for the rarity with which random spin direc- 
tions find themselves aligned. However, as this unusual 
selection bias for high spin operates only in a narrow mass 
and spin range and only for the initial detectors, we do 
not estimate this bias to an astrophysically relevant level 
of accuracy. 




FIG. 11: Detection volume for initial LIGO for BH-BH bi- 
nary populations with mi = 7712 = M/2 = IOOMq and dif- 
ferent spin orientation and magnitude distributions: aligned, 
identical spins (black, solid); aligned spins independently 
distributed uniformly between [0, a] or [— |a|,0], if a < 
(black, dotted); randomly aligned spins with identical magni- 
tude (blue, solid) ; and component spins independently drawn 
from a sphere of radius a (blue, dashed). This figure provides 
only the single-detector beaming-uncorrected detection vol- 
ume V = 47t(_D) 3 /3 assuming a single-detector SNR threshold 
p c = 8; beaming increases the detection volume by roughly 
(1.095) 3 . At low masses, far from the sharp initial LIGO 
low-frequency cutoff, all equal-mass binaries have a compara- 
ble detection volume. Bottom panel: as top panel, but now 
rri! - m 2 = 250M Q . 



VII. CONCLUSIONS 

We have provided simple phenomenological fits to the 
average detection volume within which a merging equal- 
mass BH-BH binary with M G [100, 500] Af© and ar- 
bitrary spins would produce a signal-to-noise ratio p 
greater than a fixed detection threshold, for a single 
initial LIGO interferometer operating at design sensi- 
tivity. Comparable results for advanced LIGO are de- 
scribed in the Appendix. Though we describe fits for a 
single detector, our method also applies to arbitrary net- 
works of identical detectors, adopting the universal angle- 
averaged power p combined with a suitably-modified 
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(and network-topology-dependent) beaming correction 
factor w. Our generic expression for the range is both 
surprisingly simple and accurate, involving at most 10 
(=9+1) functions to reproduce the response to generic 
spin magnitudes and orientations to cubic order in the 
component spins. Moreover, we find an even simpler 
expression involving only four mass-dependent terms 
(<^i,2,3 and ^02) also fits all our results. Though derived 
for equal mass ratio, symmetry considerations suggest 
unequal mass ratio corrections enter weakly, at higher 
order (in ex 8m x ) ■ Our spin-dependent expressions for 
p/po should therefore correctly estimate that ratio for all 
comparable-mass binaries. Finally, our method is easily 
extended to unequal mass: a i/iree-parameter symmetry- 
constrained expansion of p/p (in X± an d Sm/AI) still 
has comparatively few parameters. 

The dynamical processes that most likely produce 
merging BH-BH binaries at these masses likely guarantee 
random spin orientation [SJ. Though gravitational wave 
detectors are far more sensitive to BH-BH binaries with 
aligned spins, we conclude that BH-BH populations with 
random spin orientations will rarely provide detections 
from tightly aligned, high-mass BH-BH binaries, under 
the assumption of optimal signal processing. Lacking the 
waveforms needed to perform optimal signal processing in 
the high-mass region, however, present-day all-sky grav- 
itational wave searches conduct searches using approxi- 
mate or hierarchical methods. Further study is needed to 
assess how strongly the search methodologies themselves 
introduce bias towards aligned spin. 

Our fit suggests a surprising and astrophysically conve- 
nient conclusion: for M e [100, 5OO]M for initial LIGO 
and over M e [200, 1600]M Q for advanced LIGO, the 
population-averaged detection volume for binaries with 
random spin directions and an arbitrary spin magnitude 
distribution and a\^ < 0.8 is nearly identical (within tens 
of percent, comparable to the Poisson error in 10 detec- 
tions) to the detection volume for nonspinning binaries 
of comparable mass. On the contrary, detectors like ini- 
tial LIGO which have both a shallow optimally sensitive 
region combined with a very steep low frequency cut- 
offs will be noticably more likely to recover large aligned 
spins from a random spin population, albeit only for the 
very highest masses and spins to which the detector is 
sensitive (M > 500M o for initial LIGO). 

Conversely, our analysis implies that the detected pop- 
ulation of high-mass binaries with a given (optimally fil- 
tered) SNR p will be distributed uniformly over spin ori- 
entations, to the extent that their formation processes 
produce them. Our study therefore reaffirms the urgent 
need for models for the merger waveforms from generic 
spinning merging binaries. 

At present, our fit breaks down at very high mass, 
empirically when {p/p — 1) is of order unity. Several 
avenues of improvement could make the fit more stable: 
fitting to a logarithm hip/p , or even using similarity 
transformations to rescale p(M) to different spin geome- 
tries (e.g., changing the mass and amplitude scale using 



the ringdown frequency of the post-merger BH) . We have 
also not compared our fitting parameters with compara- 
ble coefficients for better-understood low-mass binaries 
(e.g., those undergoing simple precession), which can be 
estimated analytically and numerically. We will address 
these refinements in a future paper. 

For simplicity, our analysis is expressed in terms of 
matched filtering and a fixed SNR detection threshold. 
Whether due to incomplete signal models, inefficiencies 
in template placement, or the lack of a signal model al- 
together, realistic searches cannot identify all available 
signal power. Additionally, whether from higher mass or 
antialigned spin, waveforms of short duration are more 
easily confused with nongaussian detector noise and re- 
quire a significantly higher detection threshold. Though 
a significant technical challenge, the search- and detector- 
dependent effective detection threshold versus binary pa- 
rameters /0 C (A) that incorporates both mismatch and 
noise effects can also be tabulated. For example, for low 
mass binaries the mismatch between nonspinning search 
templates and precessing, spinning waveforms has been 
tabulated; see Brown et al. (in preparation). Combined 
with our fit to the intrinsic available detection volume, 
a "detector sensitivity" fit could efficiently communicate 
an adequately-accurate representation of the parameter- 
dependent detection volume of real high-mass searches. 

To summarize, we have provided the first phenomeno- 
logical fit to the spin-dependent detection horizon for 
generic spin magnitudes, spin orientations, and (equal) 
component masses M £ [100,800M Q ]. We have shown 
that to leading order spin effects average out of the de- 
tection volume. For example, the rate at which IMBH- 
IMBH mergers will be detected is directly proportional 
to this volume. All previous estimates for the IMBH- 
IMBH detection rate estimated this volume assuming 
minimal or occasionally aligned component spins [49Tf51] 
and therefore should be nearly unchanged even with spin 
included. 
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Appendix A: Selection biases of individual advanced 
detectors 



For pedagogical reasons, in the text we described our 
procedure estimating selection biases versus spin in the 
context of a single detector design (initial LIGO) . In this 
appendix we provide a similar discussion for advanced 
and third-generation detectors, specifically the bench- 
mark advanced LIGO |5jjj and Einstein Telescope [S3] 
designs. These instruments' sensitivity is great enough 
that, even without detections, their upper limits rule 
out otherwise astrophysically plausible progenitor mod- 
els. We particularly emphasize how these detectors su- 
perior low-frequency sensitivity leads to much more ac- 
curate fits, over a broader range of masses. 

For simplicity, despite both detectors' cosmologically 
significant range, we perform all calculations in terms of 
the luminosity distance and redshifted mass - effectively 
as if in a flat universe. Though our fits reproduce p/p 
to several percent, and though the reader can invert any 
particular line of sight to that accuracy, we have not in- 
cluded all information needed to completely reconstruct 
the selection bias versus mass and spin. At these dis- 
tances, the detection volume is no longer proportional 
to a simple cubic moment of the beampattern function; 
see for example the appendix of O'Shaughnessy et al. [I] 
for a comparable calculation at low mass. Though we 
anticipate the comoving volume swept out by the past 
detection light cone will not depend sensitively on the 
details of its truncation at large redshift 13 , we have not 
thoroughly explored the errors our neglect of these effects 
introduces. 

By adopting luminosity distance and redshifted mass 
as parameters, each result in this section is directly com- 
parable to a corresponding prediction for initial detec- 
tors. However, because advanced detectors have peak 
sensitivity at roughly 2x lower frequency, the response 
of initial detectors of mass M contains comparable wave- 
form content as and is best compared to the response of 
advanced detectors of mass 2M . For example, the low- 
mass limit for initial LIGO is roughly 1OOM for initial 
and 2OOM for advanced detectors. 



1. Average signal power versus binary spins 

Unlike the results from initial detectors, the data for 
p/po exhibits relatively weak dependence on spin for all 
masses. As a result, a fit to the numerical data per- 
forms well, both reproducing the data (Fig. |13[ ) and pro- 
ducing well-determined fitting coefficients (Fig. [l2"| ) over 
the entire range of plausible masses M 6 [200, 1600] M . 



Aside from a difference in scale, however, the fit exhibits 
properties comparable to the low- mass fit (M < 500M Q ) 
to initial LIGO p. Notably, (i) the fit is dominated by 
aligned spin coefficients, with few resolved corrections in- 
volving perpendicular spins; (ii) it depends only weakly 
on antisymmetric spin X- > (hi) it satisfies all sanity con- 
ditions, such as increasing monotonically with aligned 
spin and attaining a local extrema versus orientation 
when equal-magnitude spins aligned and antialigned; and 
finally (iv) higher harmonics I > 3 only weakly change 
the best-fit coefficient functions y a , mostly at the highest 
masses. Also as in the low-mass fit, the fitting coefficients 
are strongly correlated. Ill-constrained high-order coef- 
ficients like 23io2o are sensitive to numerical issues, such 
as extrapolation to large radius. Most differences be- 
tween the performance of initial and advanced detectors 
is directly attributable to their low-frequency response: 
advanced detectors generally lack an abrupt transition 
from peak to low-frequency sensitivity. Without a strong 
preferred frequency these detectors must have at best a 
weak bias towards recovering the largest and most tightly 
aligned spins; compare Section |IV| Furthermore, as dis- 
cussed by example in the Appendix [Bl because of the 
nearly power-law low-frequency response of the detector 
noise power spectrum, the fitting coefficients are nearly 
constant, independent of mass. 



Appendix B: Analytic estimates for expansion 
coefficients 



As seen in the text, the relative detection volume for 
different spinning binaries is largely determined by the 
angle-averaged SNR p. At both very low and very high 
mass, the leading-order dependence of the angle-averaged 
SNR on spin orientation can be calculated, particularly 
by adopting high-symmetry configurations such as spin- 
aligned binaries to better isolate relevant terms. For 
example, at very low mass the aligned spin coefficients 
such as X\ 2,3 follow from the well-known stationary- 
phase Fourier transform of the 2.5-PN-accurate inspiral 
waveform for spinning, aligned binaries. Adopting the 
restricted Newtonian amplitude at all frequencies and ex- 
panding -Jduj/dt in the denominator of the stationary- 
phase Fourier transform, the leading order change in SNR 
with spin follows approximately from 






24 J(7/3) 

III 

Si, 



(Bl) 
(B2) 



13 At large redshift the small amount of comoving volume available 
at a given redshift strongly suppresses lines of sight that reach 
back to higher redshift. 



On the other extreme, only the late-stage ringdown of 
very high-mass spinning binaries will fall into our de- 
tectors' sensitive band. Qualitatively, ringdown emission 
produces an SNR limited by the detector frequency and 
total amount of emitted power. Progenitor black hole 
spins will change the total amount of emitted energy and, 
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FIG. 12: Coefficients of fit to the ratio p/po for single design- 
sensitivity advanced LIGO detector; compare the top panel 
to Figures [8] and [8] The advanced LIGO detector design will 
be sensitive to lower characteristic frequencies [O(40 Hz) vs 
O(100 Hz) for the initial detector]. For this reason, we provide 
fits over a more limited mass range: only above 2OOM0 will 
our short waveforms start their inspiral well before they enter 
band. 



critically, the final hole ringdown frequency. Assuming 
ringdown waves are nearly monochromatic and contain 
a proportion of the total emitted energy, the high-mass 
SNR can be roughly estimated by the nonspinning am- 
plitude at that mass times an ad-hoc ringdown correction 
factor: 



P 



Po 



yflM) - M)/(M Lo - M) 

\/ Sh{frd) / Sh{frd,o) 



(B3) 



where f rd (M f ,a f ) « [1 - 0.63(1 - a / ) 3 / 10 ]/2 7 rM / is the 
ringdown frequency of the final hole of mass Mf and spin 
a/ and where Mf , a/ are known expressions of the pro- 
genitor masses and spins. If for simplicity we further 
retain only x+.z spin dependence, then we find an ex- 
pression for X\ in terms of the detector's noise power 
spectrum Sh- 



PlPo = 1 + X+,2 
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where in the last term we assume the detector's low 
frequency noise is a pure power law (Sh(f) oc /~ p ). 



FIG. 13: Expected la error T in an unconstrained fit to 
F — p/po for advanced LIGO detectors, for spin magni- 
tudes |ai | , |ci2 1 limited to below unity (blue; largest); below 
0.8 (red); and below 0.6, 0.4 and 0.2 (similar curves near bot- 
tom). Also shown (points) are the least-squares errors SF be- 
tween our fit F and the simulated NR data. Top panel shows 
standard; bottom panel shows error of restricted model. Be- 
cause of its broader bandpass, the advanced detector range is 
more accurately modeled by our expansion; compare Fig. [T] 
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FIG. 14: For an equal-mass aligned spin binary with ai = 
a2 = 0.8, estimates for p/po seen by a single advanced LIGO 
detector at D = 100 Mpc versus total binary mass M, based 
on waveforms extracted from our numerical merger simula- 
tions at coordinate extraction radii r — 40, 50, 60, 75 (blue, 
red, yellow, green) and (thick solid) linearly extrapolated in 
1/r tor-> oo. 
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FIG. 15: Coefficients of fit to the ratio p/p for single design- 
sensitivity Einstein Telescope; compare the top panel to Fig- 
ures HI and [U 



The second term is a positive-definite constant, inde- 
pendent of mass and depending predominantly on the 
physics of BH mergers; the detector influences this esti- 
mate for p only through the low-frequency noise expo- 
nent p. Though not a quantitatively precise model for 
p - the fits described in the paper are required to es- 
timate the detection volume to astrophysically-relevant 
accuracy - this simple ringdown-dominatcd signal model 
illustrates why X\ must vary relatively little between de- 
tector designs. Our estimate also identifies changes in 
d\nS h /d\nf(f = f r d(M)) as key points where X X {M) 
should change significantly. 



Appendix C: Detection volume for networks 

In the text, we estimated the detection volume to 
which a single perpendicular-arm interferometer with a 
fixed detection threshold could recover a source of un- 
known sky location. We expressed this simple measure 
of sensitivity in terms of two factors (p* , w*) . Real 
searches employ multiple detectors, providing sensitiv- 
ity to both polarizations nearly everywhere on the sky. 
Their relative orientations determine a sky position- and 
polarization- dependent sensitivity. Unfortunately, real 
searches also have detectors with different (and time- 
dependent) noise power; a network whose nodes are not 
always active; and of course beampatterns that are highly 



nonuniform. A completely realistic treatment of multidc- 
tector search sensitivity adopting realistic search strate- 
gies and thresholds is beyond the scope of this paper. 

Though necessarily incomplete, the analytic estimates 
of idealized detector performance nonetheless allow us 
to quickly understand the dominant qualitative effects 
spin might have on the detection volume, such as Fig. 
|11| Similar analytic estimates for idealized networks can 
be constructed, albeit depending on network topology. 
For coherent multidetector searches with identical detec- 
tors, their sensitivity along any direction is described in 
general by a 2x2 Hermetian matrix for each sky position 
|54j ; for brevity, we will limit attention to this case hence- 
forth. This orientation- and polarization-dependent sen- 
sitivity implies the true beampattern function w of both 
source- and sky position- orientations becomes very com- 
plex in general, particularly when multiple harmonics 
with different shapes contribute to the overall response. 
Nonetheless, there are three physically interesting ap- 
proximations where the average amplitude and beampat- 
tern are tractable: 

• Single detector (used in text): Sensitivity is char- 
acterized by an angle-averaged SNR p* [Eq. ph]; 
a beampattern function w* = p/p* of sky position 
and source orientation; and a be ampattern effective 

radius w* — (w*) [Eq. (13)], averaged over all 
orientations. 

To provide a concrete example, consider the case 
of a binary dominated by I = \m\ = 2 emission, 
viewed by a single interferometer. Since the beam- 
pattern w versus all angles is known analytically 
[e.g., Eq. (2) in [4]], the amplitude p versus orien- 
tation is completely determined by the optimally- 
oriented amplitude j0*, maX) or equivalently by the 
"horizon distance" D^ at which p max = Pc- 



p* 


P*,max 


D* 


= D h 2/5 


w* 


= 1.095 
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5 D 



/'-,„ (ci) 

(C2) 

(C3) 

Alternatively, the sensitivity of the interferom- 
eter can be characterized by source-orientation- 
averaged effective detection volume, or, to use sim- 
ilar units, the volume-averaged distance D v . The 
above expressions allow us to express familiar re- 
lationship between D v and D^ [3] in terms of the 
single-detector beampattern correction factor u)*: 



D v = D ft 2u)*/5 fa D h /(2.26) 



(C4) 



• Isotropic network, equal polarization sensitivity : 
At the other extreme from a single-detector net- 
work is the ideal case: a network with equal sensi- 
tivity to both polarizations in all directions. As- 
suming searches for this ideal network are car- 
ried out by matched filtering, its network sensi- 
tivity is characterized by an angle-averaged two- 
polarization network SNR p = p*-\/E [Eq. Q5J) ; see 
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1 VJ 10* w 

2 1.11661 1.09664 1.05065 
4 1.10614 1.08541 1.0408 
6 1.15301 1.12961 1.08489 

TABLE III: Beaming correction factors when a single I = \m\ 
mode dominates emission. In this table, w* times 1/y/E char- 
acterizes the radius of the detection volume associated with 
an all-sky search for a single source by a single interferometer; 
W (without additional factors) refers to a comparable search 
by a network with isotropic sensitivity to both polarizations; 
and w refers to a network with isotropic sensitivity to a single 
polarization. 



Cutler and Flanagan |54)]: a beampattern function 
w = p/p\ and a source-orientation-averaged mo- 
ment 



[f^fi^i/a 



„2\2/3 



(C5) 



where, because the detector has equal sensitivity to 
both polarizations, the orientation average needs to 
be conducted only over all emission directions fi. 

• Isotropic network, single polarization sensitivity: 
In between these two extremes is an isotropic net- 
work with sensitivity to only one polarization at 



each sky position. For example, the three-site 
LIGO- Virgo network is primarily sensitive to one 
polarization for each sky position. Using a sub- 
script to denote this class of detector, the the 
beampattern function can be similarly defined as a 
suitable average over emission direction fi and po- 
larization angle i\)\ 



1/3 



wo = («$) / {wl) 



(C6) 



By symmetry, the orientation-averaged signal 
power must be precisely half that of a network sen- 
sitive to two polarizations, so po — p/y2. 

To illustrate the beaming-induced differences in de- 
tection volume between these different network topolo- 



gies, in Table III we tabulate the beaming correction 
factors when a particular multipolar mode dominates 
emission (I — \m\). All are nearly unity; all depend 
weakly on I, for a given detector topology. Addition- 
ally, the beaming correction factors have a universal 
hierarchy: w > w* > w. Obviously the idealized 
two-polarization-sensitive isotropic network has the most 
symmetric beampattern, least sensitive to host orienta- 
tion and polarization content. Conversely, an isotropic 
network sensitive to one polarization is usually the most 
sensitive to unfortunately oriented (polarized) sources. 
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